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Abstract. Classical molecular dynamics enables atomistic structure simulations 
of nanoscopic systems to be made. The method is extremely powerful in solving 
the Newtonian equations of motion to predict static and dynamic properties of ex- 
tended particle systems. However, to yield macroscopically relevant and predictive 
results, suitable interatomic potentials are necessary, developed on ab-initio-based 
approximations. The fundamental requirements for performing classical molecular 
dynamics are presented as well as the relation to statistical methods and parti- 
cle mechanics, suitable integration and embedding techniques, and the analysis of 
the trajectories. The applicability of the technique is demonstrated by calculating 
quantum-dot relaxations and interaction processes at wafer-bonded interfaces. 

1 Introduction: Why Empirical Molecular Dynamics? 

Classical molecular dynamics (MD) enable atomistic structure simulations of 
nanoscopic systems and are, in principle, a simple tool to approach the many- 
particle problem. For given interatomic or intermolecular forces one has to 
integrate the Newtonian equations of motion assuming suitable boundary 
conditions for the box containing the model structure. There are at least 
two advantages of this technique. The molecular dynamics is deterministic 
and provides the complete microscopic trajectories, i.e., the full static and 
dynamic information of all particles is available, from which a large number 
of thermodynamic and mechanically relevant properties of the models can 
be calculated. Further, one can perform simulations that are macroscopically 
relevant with the present computational power of even desktop computers. 
With reasonable computational effort models of nanoscale dimension can be 
treated for several million particles and up to microseconds of real time. Thus, 
empirical MD has two main fields of application: The search for the global 
energetic minima by relaxing nanoscopic structures and the calculation of 
dynamical parameters by analyzing the lattice dynamics. 

Computer simulations are performed on models simulating the reality by 
using approximations, reduction, localization, linearization, etc., the valid- 
ity of which has to be critically evaluated for each problem considered. As 
sketched in Fig. 1, increasing the time and length scales of the models (which 
is necessary to increase the macroscopic validity and robustness of the cal- 
culations) requires an increasing number of approximations and thus leads 
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Fig. 1. Length and timescales in various modeling methods: DFT = density- 
functional theory, TB = tight-binding approximations, BOP = bond-order poten- 
tials, MD = empirical molecular dynamics, MC, CGM = Monte- Carlo/conjugate- 
gradient techniques, FEM = finite element methods, CM = continuum mechanics 



to a reduction in the ability to predict some physical properties. Density- 
functional theory (DFT) and its approximations (e.g., local-density LDA, 
generalized gradient GGA), and the different kinds of tight-binding treat- 
ments (TB) up to the bond-order potential (BOP) approximations start with 
the Born-Oppenheimer (BO) approximation, thus decoupling the ionic and 
electronic degrees of freedom. Additional approximations such as pseudopo- 
tentials, gradient corrections to the exchange-correlation potential, and in- 
complete basis sets for the single-particle states are required for DFT calcula- 
tions. The specifics of the various methods are discussed in various Chapters 
in the present book and in a number of review articles on DFT, TB, lin- 
ear scaling techniques and programs such as SIESTA and CASTEP [1-9]. In 
first-principle MD [10], e.g., using DFT or DFT-TB, the electronic system 
is treated as parameter free and the resulting Hellmann-Feynman forces are 
the glue of the ionic interactions. Such simulations are computationally too 
expensive for large systems. 

Figure 1 shows that the empirical MD closes the gap between first-princi- 
ples, macroscopic, and continuum techniques (FEM = finite element methods, 
CM = continuum mechanics). The latter neglect the underlying interatomic 
interactions but allow the description of defects, defect interactions, diffu- 
sion, growth processes, etc. Stochastic Monte Carlo techniques (MC) enable 
the further increase of the timescales, also at the first-principles level, but 
without access to the dynamics. On the other hand, static energy minimiza- 
tion (e.g., conjugate gradient methods, CGM) enable a drastic increase in 
the model size. However, this does not necessarily provide the global mini- 
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Fig. 2. Interface structures at 90° twist boundaries ((a): Pram (ra)- layer, (b): 
42m- dreidl, cf. text in Sect. 4.2), predicted by empirical MD and revealed by ab- 
nitio DFT simulations yielding semimetallic/isolating behavior as a function of the 
interface bonding state as demonstrated by the corresponding DFT band structures 
(c) and (d), respectively 

mum of the potential energy surface. The largest model dimensions (number 
of atoms) accessible to empirical MD (or MC, CGM) has approximately the 
same extension as the smallest devices in microelectronics and micromechan- 
ics. Thus, today's atomistic modeling approaches the size of actual nanoscale 
systems. 

The assumption of the existence, validity, and accuracy of known empir- 
ical interatomic potentials or force fields in analytic form always ignores the 
underlying electronic origin of the forces, i.e., the quantum structure of the 
interactions (cf. Sect. 2.4). Therefore, it is important to have better approx- 
imations, such as the bond-order potential (BOP), which is developed from 
tight-binding approximations and discussed briefly in Sect. 3.2. Other possi- 
bilities to include electronic properties of the interaction consist in continu- 
ously refitting an empirical potential during the calculation, called learning 
on the fly (see Chap. 9 and [ ]). Separated subsystems can be treated at an 
ab-initio level. Suitable handshaking methods can be designed to bridge the 
embedded subsystem with its surrounding, which is treated semiclassically 
(cf. [12-15]). However, well-constructed potentials describing sufficiently ac- 
curate physical properties also give physical insights and enable a thorough 
understanding of the underlying processes [ ]. 

Figure 2 shows a simple example demonstrating the difference in the ap- 
proximation levels. Using empirical MD the correct structural relaxation of 
special interfaces created by wafer bonding can be treated, as described later 
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in Sect. 4.2. Two configurations exist, a metastable one (Fig. 2a) and the 
global minimum structure shown in Fig. 2b, called dreidl [1 ]. The electronic 
properties demand ab-initio level simulations, done here using a smaller pe- 
riodic subunit of the interface and applying DFT techniques. It is shown 
that the higher-level approximation reproduces the structure predicted by 
the semiempirical techniques, whereas the correct energies and electronic 
properties (the band structure and the semimetallic or isolating behavior 
of the interfaces as shown in Figs. 2c, d) can only be described using the DFT 
formalism. 

A second kind of embedding problem occurs because even millions of par- 
ticles only describe a small part of reality that, even for smaller pieces of 
matter, is characterized by Avogadro's number (6 x 10 23 mol -1 ). Since iso- 
lated systems introduce strong surface effects, each model has to be embed- 
ded in suitable surroundings. For a discussion of various types of boundary 
conditions, see Sect. 3.1, Chap. 9 and [13-15]. 

The fundamental requirements of classical MD simulations, the relation to 
statistical methods and particle mechanics, suitable integration and embed- 
ding techniques and the analysis of the trajectories are presented in Sect. 2. 
The enhancements of potentials (bond-order potentials) and boundary con- 
ditions (elastic embedding) are discussed in Sect. 3. Selected application of 
semiempirical MD (relaxation of quantum dots and wafer-bonded interfaces) 
are given in Sect. 4 together with some examples from the literature. They 
strongly depend on the approximations assumed in the simulations. 

2 Empirical Molecular Dynamics: Basic Concepts 

The main steps for applying empirical molecular dynamics consist of the in- 
tegration of the basic equations (cf. Sect. 2.1) using a suitable interaction po- 
tential and embedding the model in a suitable surrounding (cf. Sect. 2.4 and 
Sect. 2.3, respectively). Textbooks of classical molecular dynamics, e.g., [17- 
23], describe the technical and numerical details, and provide a good insight 
into possible applications and the physical properties, which may be pre- 
dictable. Here, only the main ideas of empirical MD simulations, viz. the 
basic equations, methods of numerical treatment and the analysis of the tra- 
jectories are discussed. 

2.1 Newtonian Equations and Numerical Integration 

The basic equations of motion solved in empirical MD are the Newtonian 
equations for N particles (i = 1, . . . ,7V) characterized by their masses ra$, 
their coordinates r.j, and the forces acting on each particle fi. These may 
include external forces Fi and interatomic interactions fi = Y^j^i fij : 
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The assumption that a potential V exists presupposes conservative (nondissi- 
pative) interactions and needs some more general considerations, cf. Sect. 2.2. 
If one assumes pairwise central potentials, V(ri, r2, rjy) = Yli jV( r ij) 
with rij = \rj — Vi\ (the dash means that the sum is restricted to i < j), it 
follows that the virial of the external and internal forces is A4 = Yl'i j r ij fij — 
ViFi. If the system is isolated (microcanonical ensemble), the conservation 
of total energy £ = /C + V (with kinetic energy /C = rriiT 2 ) is guaranteed: 

£ = t + v = m ^ • - J2 = • ( 2 ) 

Given an initial configuration of particles and suitable boundary condi- 
tions (cf. Sect. 2.3), the differential equations can be integrated using one 
of the standard methods, e.g., Runge-Kutta techniques, predictor-corrector 
methods, Verlet or Gear algorithm. The numerical integration is equivalent 
to a Taylor expansion of the particle positions ri(t + St) at a later time in 
terms of atomic positions ri and velocities vi. The forces fi are the time 
derivatives of T{ with an increasing order at previous time steps: 

n(t + St) = n(t) + axStViit) + a 2 (St) 2 fi(t) + . . . + 0(St n ) . (3) 

In addition to good potentials and system restrictions (cf. Sect. 2.3 and 
Sect. 2.4, respectively) an efficient and stable integration procedure is re- 
quired to accurately propagate the system. The conservation of energy dur- 
ing the simulation is an important criterion. An increase in order allows 
larger time increments St to be used, if the evaluation of the higher deriv- 
atives is not too time consuming, which happens with more accurate po- 
tentials, like the BOP (Sect. 3.2). The efficiency and the accuracy of the 
integration can be controlled by choosing suitable series-expansion coeffi- 
cients dj,j = 1, 2, n — 1 and the order n of the method or by mixing the 
derivatives of Vi at different times in (3) to enhance the procedure. However, 
the increment St, of the order of fs, must be at least so small that the fastest 
particle oscillations are sufficiently sampled. 

Better or faster MD calculations may be performed using special acceler- 
ation techniques, the three most important methods being: 

1) Localization: By using the linked-cell algorithm and/or neighbor lists, 
it is assumed that for sufficiently rapidly decaying potentials (faster than 
Coulomb 1/r for r — > oo) only a small number of particles have a direct and 
significant interaction. A cutoff r c is defined and the interaction potential is 
assumed to be zero for r > r c . A transition region is fitted using splines or 
other suitable functions. Then the system is divided into cells. Their mini- 
mum dimension is given by 2r c to avoid self-interaction of the particles (min- 
imum image convention). Only the interactions of the atoms within a cell 
and its 26 neighboring cells are considered, which reduces the simulation 
time drastically from the ex N 2 behavior for all particle interactions to linear 
behavior oc 729catA/', where cat is the average number of particles per cell. 
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The problem is to find a suitable r c for a smooth transition and screening 
behavior that includes a sufficient number of next-neighbor atoms and cells. 
One also needs a suitable criterion to update the neighbor lists and reorder 
the cells whenever particles leave a cell during the propagation of the system. 

The Coulomb potential may also be screened. However, summations over 
the long-range 1/r potential, representing infinite point-charge distributions, 
are only conditionally convergent. Thus, it is better to apply the Ewald 
method (originally developed to calculate cohesive energies and Madelung 
constants [24]) and its extensions [25,26] based on successive charge neutral- 
ization by including next-neighbor shells around the origin. 

2) Parallelization: Using replicas, or dividing the structure into several 
parts, which are then distributed to different processors, thus allowing par- 
allelization [27,28]. The replica technique needs suitable criteria for dividing 
the system into small parts with minimum interaction. More importantly, 
one must bear in mind that the replicas are not independent over long times. 
A careful control of the time interval after which the communication be- 
tween the different parts is required in order to achieve the same results as 
nonparallelized simulations. This issue is the bottleneck of the technique. 

3) Time stretching: Such techniques as hyperdynamics, temperature ac- 
celeration, basin-constrained dynamics, on-the-fly Monte- Carlo, and others 
are subsumed briefly here, because their common idea consists in replacing 
the true time evolution by a shorter one increasing the potential minima, 
transition frequencies, system temperature, etc. (see [28]). 



2.2 Particle Mechanics and Nonequilibrium Systems 

In classical mechanics a system is characterized either by its Lagrangian 
£(Qi 4) — K>(Qi 4) ~ V{q) or its Hamiltonian related to the Lagrangian by 
a Legendre transform H(q,p) = ^qiPi — £ = JC + V(q), where q%^p% are 
generalized coordinates and momenta (conjugate coordinates), respectively, 
which have to be independent or unrestricted. One derives the generalized 
momenta from the derivatives of the Lagrangian pi = Jy-, whereas the deriv- 
atives of the Hamiltonian 

dU . dU 

Qi ^— , Pi -7— (4) 
dpi dqi 

reproduce Newton's law of motion (1). 

Hamilton's principle of least action enables a simple generalization of 
the mechanics of many-particle systems. The integral over the Lagrangian 
function has to be an extremum. Variational methods yields, applying the 
extremal principle: 
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The advantage of the general formulation (5) is that it allows a simple exten- 
sion to nonconservative systems by including an explicit time dependence or 
to systems with constraints, such as fixed bond lengths in subsystems, friction 
of particles, outer forces, etc. If there are holonomic constraints describing 
relations between the coordinates fi(qk) = 0, one gets the generalized ad- 
ditional coordinates aik = dfi/dqk and a\ = dfi/dt creating an additional 
term on the right-hand side of (5) in the form \ia>ik with the set of La- 
grange multipliers Xi . This formalism is the basis for the enhancement of the 
boundary conditions in Sect. 3.1 and the handshaking methods mentioned 
above. 

In addition, the Lagrangian and Hamiltonian formalisms correlate clas- 
sical dynamics to statistical thermodynamics and to quantum theory, and 
allow the evaluation of properties from the trajectories (cf. Sect. 2.5). The 
set of time-dependent coordinates q (the configuration space) and time- 
dependent momenta p (momentum space) together is called the phase space 
r = (q,p) = (qi, 92, <7n,Pi,P2, —,Pn)- Presenting f in a 6-dimensional 
hyperspace yields TV trajectories, one for each molecule, and allows the study 
of the behavior of the system using statistical methods. It is called the /i-space 
(following Ehrenfest) or the Boltzmann molecule phase space. The whole r~ 
set as one trajectory in a 6A/'-dimensional hyperspace presents the Gibbs 
phase space, with the advantage to include better interactions and restric- 
tions. However, the system must now be described as a virtual assembly for 
statistical relations. The phase space r contains the complete information on 
the microscopic state of the many-particle system. Several basic quantities 
may be derived such as the phase-space flow or the "velocity field" P = (qr, p). 
Applying the relations (4) yields the Liouville equation 




showing that the flow P behaves like an incompressible liquid, i.e., although 
the \i- trajectories are independent, their related phase-space volume is con- 
stant. 

An equivalent formulation of the Liouville statement (6) is the equa- 
tion of continuity for r, p, and the local change of the density p(q,p,i), 
similar to the Heisenberg equation for observables in quantum theory. The 
density p(q,p,t) describes the probability of finding the system within the 
region p, q and p + dp, q + dq of the phase-space volume dpdq. Because 
the systems must be somewhere in the phase space, the density can be 
normalized in the entire phase space. The integral over the whole phase 
space of the non-normalized density yields the partition function. Inte- 
gration over (N — k) particles yields the /c-particle phase-space density 
p( k \q,p,t) = J p( N \q,p,t)dp( N ~ k " > dq( N ~ k \ and the integration over the 
whole momentum space the corresponding /c-particle distribution function 
n^ k \q) = J p( N \q,p,t)dq( N ~ k ^dp. Very important for the trajectory analy- 
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sis (Sect. 2.5) are the cases k = 1,2 defining the (radial) pair-distribution 
functions g(qi,qj): 

= n^( gi , qj ). (7) 

The phase-space trajectories allow the classification of the dynamics of 
many-particle systems. They may be Poincare-recurrent (the same phase- 
space configuration occurs repeatedly), Hamiltonian (nondissipative), con- 
servative (no explicit time dependence for the Hamiltonian), or integrable 
(number of constants of motion equals the number of degrees of freedom, 
resulting in stable periodic or quasiperiodic systems). Nonintegrable systems 
may be ergodic or mixing, etc., i.e., the trajectory densely covers different 
hypersurfaces in the phase space, and allow the characterization of different 
kinds of instabilities of the system. 

2.3 Boundary Conditions and System Control 

As mentioned above, even for systems considered large in MD simulations, the 
number of atoms is small compared to real systems and therefore dominated 
by surface effects. They are caused by interactions at free surfaces or with 
the box boundaries. In order to reduce nonrealistic surface effects, periodic 
boundary conditions are applied, or the box containing the model (supercell) 
has to be enlarged so that the influence of the boundaries may be neglected. 
For further discussion and enhancements including elastic embedding, see 
Sect. 3.1. Periodic boundary conditions means that the supercell is repeated 
periodically in all space directions with identical image frames or mirror 
cells. In contrast to the case of fixed boundaries, under periodic boundary 
conditions the particles can move across the boundaries. If this happens, the 
positions r have to be replaced by r — ol where a is a translation vector to the 
image frame to which the particle is moved. Particle number, total mass, total 
energy, and momenta are conserved in periodic boundary conditions, but the 
angular momenta are changed and only an average of them is conserved. It 
should be mentioned here that periodic boundary conditions repeat defects, 
which in small supercells create high defect concentrations. Antiperiodic and 
other special boundary conditions may be chosen, where additional shifts 
of the positions and momenta along the box borders allow correction of a 
disturbed periodic continuation and a description of reflective walls at free 
surfaces to be obtained. 

According to the boundary conditions and system restrictions, the virtual 
Gibbs entities are either isolated microcanonical ensembles with constant 
volume, total energy, and particle numbers (NVE ensembles) or systems 
that exchange and interact with the environment. Closed ensembles (e.g., 
canonical NVT or isothermic-isobar NPT systems) have only an energy 
exchange with a thermostat, whereas open systems have in addition particle 
exchange with the environment (e.g., grand-canonical T/iV, where \i is the 
chemical potential). 
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Simple tools exist to equilibrate systems at a well-defined temperature. 
For example, velocity rescaling is numerically equivalent to a simple exten- 
sion of the original microcanonical system with nonholonomic constraints in 
the Lagrangian formalism and is called the Berendsen thermostat. A similar 
simple extension may be used for pressure rescaling. It is called the Berendsen 
barostat [29]. The simple velocity rescaling works by conserving the Maxwell 
distribution and yielding maximum entropy. However, such nonisolated en- 
sembles are better described using constraints in the Lagrangian similar to 
those discussed in Sect. 2.2. The extension of C in (5) by adding terms £ CO nstr 
introduces additional generalized variables having extra equations of motion 
and also additional force terms in the Newtonian equations (1). Some of the 
most important methods for £ CO nstr are given below without further com- 
ments: 

- Nose-bath and generalized Nose-Hoover thermostat [30-32]: £ C onstr = 
Ms 2 / 2 — rifTln(s) with the fictive or the whole mass M = J] m ^ the 
degrees of freedom rif = 3N + 1, and the new generalized variable s 
playing the role of an entropy. 

- Andersen isothermic-isobaric system control (NPH ensemble) £ CO nstr = 
V 2 /2 J r PV introduces volume and pressure as generalized variables [33]. 

- Generalized stresses according to Parrinello and Rahman [34] NTLa- 
ensembles: the Lagrangian is extended by £ C onstr = — l/2TrXJG with a 
generalized symmetric tensor U and the metric tensor G of the crys- 
tal structure. A further specialization, e.g., for constant strain rates 

Gyy &zz -ensemble [35], is in principle a mixture of the isothermic- 
isobaric system with the generalized stress constraint. 

- Brownian fluctuations, transport and flow processes, density and other 
gradients, Langevin damping ex v(t), etc. demand nonequilibrium MD 
[36, 37], which is mostly done by including the perturbation as a suit- 
able virtual field J-(p, q, t) into the equations of motion via the Lagrange 
formalism. 

2.4 Many-Body Empirical Potentials and Force Fields 

Empirical potentials and force fields exist with a wide variety of forms, and 
also different classification schemes are used according to their structure, ap- 
plicability, or physical meaning. It makes no sense to describe a large number 
of potentials or many details in the present review, therefore only some of the 
existing and most used potentials are briefly listed below (a good overview 
can be found in [38-48]. A classification by Finnis [2] (and many references 
therein) describes recent developments and discusses in an excellent way the 
justification of the potentials by first-principle approximations, which is im- 
portant for the physical reliability and the insight into the electronic structure 
of potentials. Interatomic forces are accurate only if the influence of the lo- 
cal environment according to the electronic structure is included. The most 
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important property of a potential is transferability, that is the applicability 
of the potential to varied bonding environments. An important additional 
criterion is its ability to predict a wide range of properties without refitting 
the parameters. The number of fit parameters decreases as the sophistication 
of the force fields increases. According to [ ] one has: 

Pair potentials - Valid for s-p bonded metals and mostly approximated 
by a sum of pairwise potentials. They may be derived as the response 
to a perturbation in jellium, which can be visualized in the pseudoatom 
picture as an ionic core and a screening cloud of electrons. 
Ionic potentials - To use Coulomb interactions directly for ionic struc- 
tures, the problem of screening (Ewald summation, Madelung constant) 
has to be considered as mentioned above. The interactions were originally 
described by Born and others as the rigid-ion approximation. Starting 
from the Hohenberg-Kohn-Sham formulation of the DFT-LDA, shell or 
deformable ion models may be developed beyond the rigid-ion approx- 
imation, where the additional shell terms [49] look like electron-density 
differences in noble gases. 

Tight-binding models - Different derivations and approximations for TB- 
related potentials exist and are nowadays applicable to semiconductors, 
transition metals, alloys and ionic systems. The analytic bond-order po- 
tential is such an approximation and will be discussed in more detail in 
Sect. 3.2. 

Hybrid schemes - Combinations of pair potentials with TB approxima- 
tions are known as generalized pseudopotentials, effective media theories 
(EMT [ ]), environmental-dependent ionic potentials (EDIP [ ]), or em- 
bedded atom models (EAM [52]); for details see [2]. 

From the empirical point of view the simplest form of potentials may 
be considered to be a Taylor-series expansion of the potential energy with 
respect to 2-, 3-,. . ., n-body atomic interactions. Pair potentials have a short- 
range repulsive part, and a long-range attractive part, e.g., of the Morse or 
Lennard-Jones (LJ) type, mostly "12-6", i.e., ar~ 12 — br~ 6 . LJ potentials 
are successfully applied to noble gases, biological calculations, or to model 
long-range van der Waals interactions (e.g. [53]). For quasicrystals [54] an 
LJ potential was constructed with two minima in a golden number distance 
relation. However, simple pair potentials are restricted in their validity to 
very simple structures or to small deviations from the equilibrium. Therefore 
many-body interactions are added and fitted for special purposes, e.g., the 
MD of molecules and molecule interactions [55,56]. Separable 3-body interac- 
tions are widely used: Stillinger- Weber [57] (SW), Biswas-Hamann [58], and 
Takai-Halicioglu-Tiller [59] potentials. The SW is perhaps the best-known 
3-body-type potential. It includes anharmonic effects necessary to reproduce 
the thermal lattice expansion of Si and Ge [60] . Hybrid force fields are some- 
times used to include the interaction of different types of atoms, such as Born- 
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Mayer-Huggins, Rahmann-Stillinger-Lemberg terms and others applicable 
to silicate glasses and interdiffusing metal ions or water molecules [61-65]. 

As mentioned above, other force fields are developed from first-principle 
approximations that combine sufficient simplicity with high rigor. They are 
not based upon an expansion involving TV-body interactions (cluster poten- 
tials). The more or less empirical forms of TB potentials and effective medium 
force fields are the modified embedded atom model (MEAM, [ ] for cubic 
structures and references therein for other structures), the Finnis-Sinclair 
(FS) [ ] and the Tersoff-type (TS) potentials [68-70]. The TS potential is 
an empirical bond-order potential with the functional form: 



The bonds are weighted by the bond order bij = F(rik,rjk,'Yijk) including 
all next neighbors k ^ i,j, which gives the attractive interaction the form 
of an embedded many-body term. The different parameterizations (TI, Til, 
Till) of the Tersoff potential have been intensively tested. Other parameteri- 
zations exist [71]. They involve other environmental functions and first-prin- 
ciple derivations [72,73], as well as extensions to include further interactions, 
H in Si [ ], C, Ge [75], C-Si-H [76,77], AlAs, GaAs, In As, etc. [78]. Finally, 
a refitted MEAM potential with SW terms is available for Si [ ]. Multi- 
pole expansions replaced by spherical harmonics [80] are an alternative to 
TS potentials. 

A comparative study of empirical potentials shows advantages and disad- 
vantages in the range of validity, physical transparency, fitting and accuracy 
as well as applicability [81]. Restrictions exist for all empirical potential types, 
even if special environmental dependencies are constructed to enhance the 
elastic properties near defects. In addition, not all potentials are applicable 
to long-range interactions, and the electronic structure and the nature of the 
covalent bonds can only be described indirectly. Thus, it is of importance to 
find physically motivated semiempirical potentials, as mentioned above and 
discussed in Sect. 3.2 for TB-based analytic BOPs. The parameters of the 
empirical force field have to be fitted to experimental data or first-princi- 
ple calculations. First, the cohesive energy, lattice parameter and stability of 
the crystal structures have to be tested or fitted. The bulk modulus, elastic 
constants, and phonon spectra are very important properties for the fit. The 
following section describes some of the quantities that may be used for the 
fit or to be evaluated from the MD simulations if not fitted. Very important 
details concerning point defects and defect clusters - necessary to get the 
higher-order interaction terms - are given by the energy and structure of 
such defects. The data may be given by DFT or TB dynamics or geometry 
optimizations, as, e.g., [82-85]. 



V(rij) = ae 



-Xr, 




(8) 
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2.5 Determination of Properties 

Static properties of systems simulated by empirical MD can be directly cal- 
culated from the radial- or pair-distribution function (7). Dynamic properties 
follow from the trajectories using averages or correlations: 



which in principle correspond to time averages. However, they are sampled at 
discrete points, so that it is necessary to choose suitable sampling procedures 
to reduce the effects of the finite size of the system, stochastic deviations, 
and large MD runs. 

Two basic relations are central for the analysis of the properties. The 
ergodic hypothesis states that the ensemble average (A) e is equal to the 
time average (A) tl which relates the averages to the measurement of a sin- 
gle equilibrium system. The Green-Kubo formula lim^oo ^ A( ^~_^ to ^ - = 

dr((A(r)A(to)) e ) relates mean square deviations with time correlations. 
The diffusion coefficient D = lim^oo — r j(0)] 2 ) (Einstein rela- 

tion) is equivalent to the velocity autocorrelation function, which is a special 
form of the Green-Kubo formula D = / °°(^ v j(t) ' v j(®))- Similarly, one 
uses the crosscorrelation of different stress components to obtain shear vis- 
cosities. Other transport coefficients may be derived analogously. 

In thermodynamic equilibrium the kinetic energy JC per degree of freedom 
is determined by the equipartition theorem JC = Ei m ^ 2 /2) = 3Nk^T/2 
(&b = Boltzmann constant) which yields a measure of the system temper- 
ature T. The strain tensor a and the pressure P are obtained from the 
generalized virial theorem a M = V^Ej v jk • v jt + J2ij r ijk • fiji]- Thus, 
the pressure is the canonical expectation value P = 1/3F[2/C — M) of 
the total virial A4. Alternatively, one can use the pair distribution in the 
form P = pk^T — p 2 J °° g(r) ^47rr 3 dr, which may be useful for correcting 
sampling errors. 

Using the densities p(r,p) as defined above and Z = J p(r,p)dr as nor- 
malization, statistical mechanics deals with ensemble averages, which in gen- 
eral are written as 



All thermodynamic functions may be derived from the partition func- 
tion Z, for example the free energy F(T,V,N) = -k B T\n[Q(T, V, N)]. In 
addition, one can obtain all the thermodynamic response coefficients. With 
the internal energy U, one computes the isochoric heat capacity Cy = 
{dU/dT)jsfy and thus a measure for the quality of the temperature equi- 
libration ((T) 2 - (T 2 ))/(T 2 ) = 3/2iV(l - 3kN)/2c v . From the volume V 
follows the isothermal compression xt — — 1 /V \dV / dP) n ,t , etc. One has to 
choose according to the different ensembles: 




(9) 




(10) 
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- Microcanonical: Pnve = S[H(r,p) — E] conserving entropy S(E,N, V). 

- Canonical: Pnvt — e^ H ^ kuT ^ oc e v ^ r ^ kBT conserving free energy F(T, 



- Isothermic-isobaric: Pnpt — ^ H TS ^ k ^ T conserving Gibbs free energy 



- Grand canonical: pTfiV = ^ H ^ N ^ k ^ T conserving Massieu function 
J(T, /i, V) (Legendre transformation in entropy representation). 

Finally, it should be mentioned that the Fourier transform of pair dis- 
tributions is connected to the scattering functions in X-ray, neutron and 
electron diffraction. MD-relaxed structure models allow the simulation of the 
transmission electron microscope (TEM) or high-resolution electron micro- 
scope (HREM) image contrast and therefore make the contrast analysis more 
quantitative. For this purpose, snapshots of the atomic configurations are cut 
into thin slices, which are folded with atomic scattering amplitudes and each 
other to describe the electron scattering (multislice formulation of the dynam- 
ical scattering theory), cf. the applications in Sect. 4 using this technique to 
interpret HREM investigations of quantum dots and bonded interfaces. 

3 Extensions of the Empirical Molecular Dynamics 

Coupling of length and timescales in empirical MD means bridging the first- 
principles particle interactions and the box environment (1). It can be done 
either using embedding and handshaking or by a separate treatment and 
a transfer parameter between the subsystems. MD simulations of the crack 
propagation [86] and the analysis of submicrometer MEMS [ ] are success- 
fully applications of the FEM coupling between MD and an environmental 
continuum. In Sect. 3.1 enhanced boundary conditions for MD are discussed: 
where the coupling between MD and an elastic continuum is a handshaking 
method based on an extended Lagrangian [88, 89]. The main steps in the 
development of an analytic TB-based BOP [90] are sketched in Sect. 3.2 as 
an example of using enhanced potentials. 

3.1 Modified Boundary Conditions: Elastic Embedding 

Elastic continua may be coupled to MD when the potential energy of an 
infinite crystal with a defect as shown in Fig. 3a is approximated in the outer 
region II by generalized coordinates a& [89] : 



In the defect region I, characterized by large strains, the positions of atoms 
Vi (i = 1, . . . , N) are treated by empirical MD. The atomic positions rj (j > 
N) in the outer regions II and III result from the linear theory of elasticity 



V,N) 



G(T,P,N) 



E({r i },{r j }) = E({r i },{a k }). 



(11) 



rj = Rj + u^>(Rj) + u(Rj, {a k }) 



(12) 
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WW 

Fig. 3. Dislocation geometry ((a), I = MD region, II = elastic, III = overlap, 
cf. text) to apply elastic boundary conditions for a dipole of 60° dislocations, and 
snapshots during MD annealing: 500 K (b), 600 K (c), OK (d) 




where the fields u°(R) and u(R) describe the displacements of atoms from 
their positions Rj in the perfect crystal, and satisfy the equilibrium equations 
of a continuous elastic medium with defects; u°(R) is the static displacement 
field of the defect and independent of the atomic behavior in I; u(R) is 
related to the atomic shifts in region I and can be represented as a multipole 
expansion: 



u(R,{a k }) = ^a k uW(R) 



(13) 



k=i 



over homogeneous eigensolutions U^ k \R) [ ..]. The are rapidly decreas- 
ing with R, thus the sum in (13) is truncated to a finite number K of terms. 

The equilibrium positions of atoms in the entire crystal are obtained from 
the minimization of the potential energy given in (11) with respect to r{ 
and ajfe. This is equivalent to a dynamic formulation based on the extended La- 
grangian as discussed above (Sect. 2.2 and Sect. 2.3) with the extension (11): 



N . 2 

m i r j 



K • 2 

jJ>Q>k 



E({r z },{a k }). 



(14) 



k=i 



Here, rrii are the atomic masses and \i is a parameter playing the role of 
mass for the generalized coordinates a k . If \i is properly chosen, the phonons 
are smooth from I to II and the outer regions oscillate slower than the MD 
subsystem as demonstrated in the snapshots Figs. 3b-d. The corresponding 
system defining the forces reads 



dE({n},{a k }) 



(15) 



dE({ ri },{a k }) 
da k 



j>N 



dE({ ri },{rj}) drj 



drj 



da k 



^FjU^iRj). 

jen 

(16) 
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These must vanish in the equilibrium Fi = F k = 0. The sum in (16) expands 
only over the atoms in region II surrounding region I, since in region III the 
linear elasticity alone is sufficient and the forces on these atoms are nearly 
vanishing. However, region II must be completely embedded in region III. The 
forces on atoms in region II must be derived from the interatomic potential, 
therefore the size of region II must be extended beyond the potential cutoff r c . 
The Lagrangian (14) results in the equations of motion 

mr i = F i , (i = l,...,N), n'd k = F k , (k = l,...,K), (17) 

thus extending Newton's equations by equivalent ones for the generalized 
coordinates. 

3.2 Tight-Binding-Based Analytic Bond-Order Potentials 

As discussed before, the use of TB methods allows much larger models than 
accessible to DFT. TB formulations exist in many forms (e.g., [2, 6-8, 91- 
93]). The application of an analytical BOP, developed mainly by Pettifor 
and coworkers [90, 94-98], may further enhance empirical MD by providing 
better justified force fields while allowing faster simulations than numerical 
TB-MD. The approximations to develop analytic TB potentials from DFT 
may be summarized by the following steps (cf. also [2,99]): 1. Construct the 
TB matrix elements by Slater-Koster two-center integrals including s- and 
p-orbitals, 2. transform the matrix to the bond representation, 3. replace the 
diagonalization by Lanczos recursion, 4. get the momenta of the density of 
states from the continued fraction representation of the Green's function up 
to order n for an analytic BOP-n potential. The basic ideas sketched with 
a few more details are as follows. 

The cohesive energy of a solid in the TB formulation can be written in 
terms of the pairwise repulsion f7 rep of the atomic cores and the energy due 
to the formation of electronic bonds 

U C oh — U Te p -\~ t^band ^atoms 

= ^£V(ry ) + 2 £ e<»> - £ N^ m e ta , (18) 

i,j n(occ) iat 

where the electronic energy e of the free atoms has to be subtracted from the 
energy of the electrons on the molecular orbitals e. Replacing with the 
eigenstates of the TB-Hamiltonian for orbital a at atom i, 

e (») = e ^(n\n) = (n\e™\n) = (n\H\n) = £ C% ] H ia ^C% ] , (19) 

ia,jP 

the electronic contributions to the cohesive energy L^band — ^atom = ^bond + 
^prom can be rearranged and separated in the diagonal and off diagonal parts: 

£4>ond = 2 ^ Pjf3,iaHia,j{3, U prom = ^ [%pia,ia ~ ^i^° m ] ^iot • (20) 
ia,j{3 iot 
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The third contribution, the promotion energy U prom compares the occu- 
pancy of the atomic orbitals of the free atoms with the occupation of the 
corresponding molecular orbitals [100]. The bond energy Uh on d describes the 
energy connected with the exchange or hopping of electrons between arbitrary 
pairs of atomic neighbors the factor 2 is due to the spin degener- 

acy. The transition matrix elements of the Hamiltonian are hopping energies, 
and their transition probability is given by the corresponding element of the 
density matrix. The contribution for one bond between the atoms i and j 
is thus characterized by the part of the density called the bond order Ojp^a 
that may be expressed as a trace: 

^b'id = E "--i '■■>■■ = Tr{H9) . (21) 

a,f3 

Besides the bond energy, the force exerted on any atom i must be given 
analytically and therefore one needs the gradient of the potential energy 
in (18) at the position of the atom i: 

dU coh dU hond 1 ^ d(j){rij) 



dri dri 2 ^— ' dri 

j 

This expression includes electronic bonds and ionic pairwise repulsions from 
all atoms of the system. The general form is still expensive to cope with for 
the simulation of mesoscopic systems. O (TV) scaling behavior is provided if 
the cohesive energy of the system is approximated by the Tersoff potential 
in (8), where bije~^ Vi3 is replaced by U^ nd an( ^ a suitable cutoff with resulting 
balanced interatomic forces is added. 

To find an efficient way to obtain the bond energy in a manner that scales 
linearly with the system size, too, the derivative dU Q° nd of the bond energy 
in (22) is replaced assuming a stationary electron density p in the electronic 
ground state and applying the Hellmann-Feynman theorem. The forces can 
now be obtained without calculating the derivatives of the electronic states 
and leads to the Hellmann-Feynman force [101, 102]: 

- f "=EV^. (23) 

ia>,j(3 

Both the elements of the density matrix and the hopping elements of the 
Hamiltonian are functions of the relative orientation and separation of the 
bonding orbitals and have been calculated by Slater and Koster [ ] as- 
suming a linear combination of atomic orbitals (LCAO). To introduce the 
Slater-Koster matrix elements, usually denoted by ssa, spa,pp7r, etc., corre- 
sponding to the contributing orbitals and its interaction type, the TB Hamil- 
tonian operator is transformed to a tridiagonal form. This can be done using 
the Lanzcos transformation [104]: 

(U m H)U n (24) 

H\U n ) =a n \U n ) + b n \U n - 1 )+b n+1 \U n+1 ), (25) 
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where the orthonormal basis at the higher level n is recursively developed 
from \Uq). The coefficients a n , b n are the elements of the continuous fraction 
of the Green's function below. 

The offdiagonal elements of the density matrix are related to the Green's 
function [105], 

1 f Ef 

piajp = - - lim 3? / dE G iaJ p(Z) , (26) 

where the complex variable Z = E + irj is the real energy E with a positive, 
imaginary infinitesimal to perform the integration (theorem of residues). This 
intersite Green's function can be connected to the site-diagonal Green's func- 
tion [105] via 

Giajp(Z) = d^ 00 ^ ^ + GQ (Z)S ii jSa 1 /3 , (27) 

and the latter can be evaluated recursively [106] using the coefficients of the 
Lanczos recursion algorithm as mentioned above: 

GUZ) = 1 —[7^ • (28) 

7-„A \2) 



The bond order can now be expressed in terms of the derivatives of the 
recursion coefficients a£ and b£. 



®ia>,j(3 — 2 



f)n A f)h A 

V 'V A n 4- ? V w A n 

X0n,n0 ~ a + Z ^0(n-l),n0 a, . 

n=0 n=l 



(29) 



with the response function Xom,no(Ef) = - Hm^^o ^ J Ef dEGQ m (Z)G^ (Z) 
and the elements Go n = G n o following from the system of equations 
(Z - a n )G nm (Z) - b n G n - ly7n (Z) - b n+1 G n+ljrn (Z) S ny7n . The more re- 
cursion coefficients included in (29), the more accurately the bond order will 
be approximated. The recursion coefficients are related to the moments of the 
local density of states (LDOS) [105] and the site-diagonal Green's function 
of (26) and (27) relates to the LDOS itself. Therefore, the recursive solution 
of (28) gives an approximation to LDOS in terms of its moments [107] 



(n) 



J E n n ia (E)dE = ^2 H ia j l(3l H jlf3l j 2( 3 2 .. . H jn _ l(3n _ liia (30) 



and may be interpreted as self-returning (closed) loops of hops of length n 
for electrons over neighboring atoms. The local atomic environment defines 
the LDOS via the moments (30), which in turn is used to calculate the bond 
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order in (21) and the (local) atomic force (23). The remaining free parame- 
ters in the analytic form (8) with U^ nd instead of bij may be fitted in the 
usual way. This is still a hard task, because the bond and the promotion 
energy involve different parameters as the repulsive part and cutoff parame- 
ters and screening functions for all terms have to be included. Besides an 
accurate fit, the BOP requires well-parameterized TB matrix elements or pa- 
rameter optimizing, and the problem of transferability [99, 108-110] has to 
be considered separately. For BOP of order n = 2 the bond-order term in the 
TS-representation reads = (—ssciij + pp(Jij)0ia,ja — ^PP^ijOinj-K and the 
numerical behavior of BOP 2 and TS are approximately equivalent. The de- 
tails for higher-order BOP are given in the papers of Pettifor's group [2, 99]. 
The b^ terms of the analytic BOP4 involve complex angular dependencies, 
partially beyond those neglected in Pettifor's formalism. For structures with 
defects as well as the wafer bonding of diamond surfaces where 7r-bonds can- 
not be neglected, BOP can be found in [111-114] and will soon be published 
elsewhere with detailed derivations. 



4 Applications 

It is impossible to review the rapidly growing number of successful applica- 
tions of empirical MD in materials research. A few representative examples 
may give an impression of the wide range of problems considered. Isolated 
point defects are mainly simulated to check the quality of the fit of the poten- 
tial parameters. Surface reconstructions, adatom and absorption phenomena, 
growth processes, and especially extended defect structures and interactions 
can be investigated in detail. The analysis of dislocation core-structures [115], 
the use of core-structure data for studying dislocation kink motion [116] or 
the dislocation motion during nanoindentation [117] are examples. Interface 
investigations have a long tradition, as illustrated in the standard book for 
grain-boundary structure [ ] and growth [119]. Heterophase interfaces, e.g., 
using a Khor-de-Sama potential for Al and TS for SiC [12 ] to simulate 
Al/SiC interfaces, demand special attention to the correct description of the 
misfit [121] to get good interface energies. Simulations with the Tersoff po- 
tential and its modifications yield the correct diameter of the critical nuclei 
for the growth of Ge nanocrystals in an amorphous matrix [122], allow the 
study of growth, strains, and stability of Si- and C-nanotubes [123, 124], 
and SiC surface reconstructions to propose SiC/Si interface structures [125]. 
A review of atomistic simulations of diffusion and growth on and in semi- 
conductors [126] demonstrates the applicability of SW and TS potentials in 
comparison with data from TB and DFT calculations. 

In the following, two examples of our work are discussed. Empirical MD 
simulations of first, high-resolution electron microscopy (HREM) image con- 
trasts in quantum dots (Sect. 4.1) and second, the physical processes at in- 
terfaces during wafer bonding (Sect. 4.2). 
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4.1 Quantum Dots: Relaxation, Reordering, and Stability 

A quantum dot (QD) is a nanometer-scaled island or region of suitable ma- 
terial free-standing on or embedded in semiconductor or other matrices. The 
possibility to arrange QDs into complex arrays implies many opportunities for 
scientific investigations and technological applications. However, depending 
on the growth techniques applied (mainly MBE and MOCVD), the islands 
differ in size, shape, chemical composition and lattice strain, which strongly 
influences the confinement of electrons in nanometer-scaled QDs. The shape, 
size and strain field of single QDs, as well as quality, density, and homogeneity 
of equisized and equishaped dot arrangements determine the optical proper- 
ties, the emission and absorption of light, the lasing efficiency, and other 
optoelectronic device properties [127, 128]. A critical minimum QD size is re- 
quired to confine at least one electron/exciton in the dot. A critical maximum 
QD size is related to the separation of the energy levels for thermally induced 
decoupling. Uniformity of the QD size is necessary to ensure coupling of states 
between QDs. The localization of states and their stability depend further on 
composition and strain of the QDs. The strain relaxation at facet edges and 
between the islands is the driving force behind self-organization and lateral 
arrangement, vertical stacking on top or between buried dots, or preordering 
by surface structuring. In addition, an extension of the emission range to- 
wards longer wavelengths needs a better understanding and handling of the 
controlled growth via lattice mismatched heterostructures or self-assembling 
phenomena (see, e.g., [129]). 

A wide variety of imaging methods are used to investigate growth, self- 
assembly, and physical properties of quantum dots. Among these the cross- 
sectional HREM and the plan-view TEM imaging techniques are suitable 
methods to characterize directly shape, size, and strain field [ ]. But the 
HREM and TEM techniques images are difficult to interpret phenomeno- 
logically, especially when separating shape and strain effects. Modeling is 
essential to uniquely determine the features and provide contrast rules. 

In Fig. 4, the atomic structure of an InGaAs-QD in a GaAs matrix (for 
other systems cf. [130, 131]) is prescribed by geometric models and relaxed 
by MD simulations. Very different dot shapes have been proposed and the- 
oretically investigated: lens-shaped dots, conical islands, volcano-type QDs, 
and pyramids with different side facets of type {Oil}, {111}, {112}, {113}, 
{136}, and both {011} + {111} mixed, etc. Some of these and a spherical cap 
are schematically presented in Fig. 4a; in simulations one or two monolay- 
ers (ML) thick wetting layers are included, too. The most important differ- 
ence between the various structures is the varying step structure of the facets 
due to their different inclination. There are at least two reasons to inves- 
tigate these configurations [130]. First, small embedded precipitates always 
have facets; a transition between dome-like structure and pyramids due to 
changes in spacer distance, change the number and arrangement of the facets, 
and thus strain and electronic properties. Second, for highly faceted struc- 
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m m n 121 *■ (c) non relaxed relaxed (d) 



Fig. 4. Structure modeling and image simulation of different pyramidal-shaped 
quantum-dot configurations: (a) different faceting, truncation, and wetting of pyra- 
midal start models (matrix removed, models related to (OOl)-base planes), (b) re- 
laxed complete model of a {011} pyramid (base length about 6nm, 10 x 10 x 10- 
supercell length 10 nm), (c) energy relaxation of a {011} quantum dot (poten- 
tial E po t and total Etot energy versus time steps), (d) cross-sectional HREM 
and bright-field diffraction contrast simulated for model (b) before and after 
relaxation assuming a standard 400 kV microscope at the Scherzer focus (i.e., 
A -40 nm, C s = 1mm, a 1.2 nm" 1 , t = 9 nm, 5 = 8nm"\ cf. [130]) 
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tures the continuum elasticity is practically inapplicable and finite element 
calculations must be done in 3-D instead of 2-D. However, the technique of 
ab-initio MC provides a more accurate prediction of the QD shapes than em- 
pirical MD [132]. The embedding of one perfect {011} pyramid in a matrix is 
demonstrated in Fig. 4b after prerelaxation. Figure 4c shows typical anneal- 
ing behavior during empirical MD calculation, characterized by the potential 
E pot and the total energy E tot per atom. The energy difference E^ot — Epot 
is equal to the mean kinetic energy and is thus directly related to the tem- 
perature of the system. After the prerelaxation of 5 ps at K, an annealing 
cycle follows, 60 ps stepwise heat up to 600 K and cool down to OK, equi- 
librating the system at each heating step. The example here demonstrates 
a short cycle; most of the embedded QDs are relaxed at each T-step for at 
least 10 000 timesteps of 0.25 fs, followed by annealing up to about 900 K 
(the temperature is not well defined with empirical potentials but is below 
the melting temperature). Whereas the structure in Fig. 4b is less strained, 
highly strained configurations occur due to the self-interaction of the QD in 
small supercells that correspond to a stacked sequence with very small dot 
distances. The extension of the supercell chosen in the simulations depends 
on the extension of the QD to be investigated, as well as on the overlap of 
the strain fields at the borders, to avoid self-interactions. Whereas Fig. 4b 
shows a relatively small supercell, investigations are made for supercells of up 
to 89 x 89 x 89 unit cells with a base length of the QD of 9 nm containing sev- 
eral million atoms. For the image simulations, subregions of the supercells are 
used, sliced into one atomic layer and applying the multislice image simulation 
technique. By comparing imaging for structures before and after relaxation, 
Fig. 4d demonstrates the enormous influence of the relaxations on the image 
contrast in cross-sectional HREM and TEM [130, 131]. In summarizing some 
of the results one can state that MD calculations with SW (CdZnSe) and 
Til (InGaAs, Ge) allow us to obtain well-relaxed structures, in contrast to 
the static relaxations performed with the Keating potential [133, 134]. With 
this insight into the atomic processes of rearranging and straining QDs at 
an atomic level the growth conditions for quantum dots may be enhanced as 
a first step to tailoring their properties. 

4.2 Bonded Interfaces: Tailoring Electronic 
or Mechanical Properties? 

Wafer bonding, i.e., the creation of interfaces by joining two wafer surfaces, 
has become an attractive method for many practical applications in micro- 
electronics, micromechanics or optoelectronics [135]. The macroscopic prop- 
erties of bonded materials are mainly determined by the atomic processes 
at the interfaces (clean and polished hydrophobic or hydrophilic surfaces, as 
schematically shown in Fig. 5) during the transition from adhesion to chemi- 
cal bonding. For this, elevated temperatures or external forces are required, as 
could be revealed by MD simulations of hydrogen-passivated interfaces [136] 
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Fig. 5. Wafer bonding at Si [100] surfaces: simulations of surface rearrangement 
for perfect alignment and slow heating. The figure describes the evolution of the 
atomic relaxation leading to bonded wafers 

or of silica bonding [ 5]. Thus, describing atomic processes is of increasing 
interest to support experimental investigations or to predict bonding behav- 
ior. Already, slightly rescaled SW potentials predict the bond behavior via 
bond breaking and dimer reconfiguration as shown in the snapshots explain- 
ing the possibility of covalent bonding at room temperature for very clean 
hydrophobic surfaces under UHV conditions [137]. 

Whereas the bonding of two perfectly aligned, identical wafers gives a sin- 
gle, perfectly bonded wafer without defects, miscut of the wafer results in 
steps on the wafer surfaces and edge dislocations at the bonded interfaces 
are created. In Fig. 6 the red color describes the potential energy above the 
ground state during wafer bonding over two-atomic steps. The upper ter- 
races behave like perfect surfaces and the dimerized starting configuration 
of Fig. 6a is rearranged and new bonds are created. The energy gained is 
dissipated and for slow heat transfer the avalanche effect leads to bonding 
of the lower terraces, too. Two 60° dislocations remain and, depending on 
the rigid shift of the start configuration, an additional row of vacancies [137]. 
The dislocations may split as simulated and observed in experiments [138]. 
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Fig. 6. MD simulation of wafer bonding over surface steps. Note the dislocation 
relaxation for the initial configuration (left), a snapshot during annealing up to 
900 K, after 12.5 ps (middle) and the final relaxed state (right) 

Bonding wafers with rotational twist leads to a network of screw disloca- 
tions at the interface. A special situation is the 90° twist, which always occurs 
between monoatomic steps. A Stillinger- Weber potential applied to a 90°- 
twist bonded wafer pair [16] yields a metastable fivefold-coordinated interface 
with a mirror symmetry normal to the interface characterized by a Pmm(m)- 
layer group (cf. Fig. 2a). Using the Tersoff or BOP-like potentials [113] and 
metastable or well-prepared starting configurations allows further structure 
relaxation and energy minimization. Figure 2b shows this relaxed configu- 
ration, which is (2 x 2) reconstructed and consists of structural units with 
a 42m- (D2d) point group symmetry, called the 42m- dreidl. It should be em- 
phasized that the dreidl structure is found to be the minimum energy config- 
uration also in DFT-LDA simulations [16]. However, the energies differ from 
those given in [139, 140]. Much more important is the modification of the 
band structure due to the different interface relaxation that may enable en- 
gineering of the electronic properties: whereas the metastable configuration 
(cf. Fig. 2c) yields semimetallic behavior, the dreidl structure (cf. Fig. 2d) 
yields a larger bandgap than in perfect lattices. The dreidl interface structure 
and its band-structure modification is very similar to the essential building 
blocks proposed by Chadi [141] for group-IV materials. They describe geom- 
etry and properties of the transformation of Si and Ge under pressure and 
the special a//o-phases as a new class of crystalline structures. 

A small misalignment of the wafers during wafer bonding yields bonded 
interfaces with twist rotation resulting in a checkerboard-like interface struc- 
ture [142, 143]. Figure 7 shows some of the resulting minimum structures 
gained for higher annealing temperatures and different twist rotation angles. 
Before the bonding process takes place, the superposition of the two wafers 
looks like a Moire pattern in the projection normal to the interface. After 
bonding and sufficient relaxation under slow heat- transfer conditions, almost 
all atoms have a bulk-like environment separated by misfit screw disloca- 
tions, which may have many kinks. The screw dislocation- network of the 
bonded wafer has a period half that of the Moire pattern. One finds more 
imperfectly bonded regions around the screw dislocations for smaller twist 
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Fig. 7. MD relaxations of bonding rotationally twisted wafers ([001] and [110] 
views) with different angles: (a) 2.8°, 22 nm box, orthogonal dimers, [001] view; 
(b) 6.7°, 9.2 nm box, orthogonal dimers, [001] and [110] views; (c) as (a) and (d) as 
(b) with parallel dimers 

angles, whereas bonding at higher angles results in more or less widely spread 
strained interface regions. In Fig. 7 bonding with orthogonal start configu- 
rations (a,b) is compared to those with parallel dimer start configurations 
(c,d). Thus the bonding of Figs. 7a and b may be considered as bonding with 
an additional 90° twist rotation. Clearly the periodicity of the defect region 
is twice those of Figs. 7b and c smoothing out the interface, but creating ad- 
ditional shear strains. Irrespective of the chosen twist angles and box dimen- 
sions all final structures yield bond energies of approximately 4.5eV/atom 
at OK, however, varying slightly with the twist angle. A maximum occurs 
between 4° and 6° twist related to a change of the bonding behavior itself. 
The higher the annealing temperature the better the screw formation. In con- 
clusion, simulation of the atomic processes at wafer-bonded interfaces offers 
not only the tailoring of electronic properties, it is also an important step in 
understanding how to control the creation of special interface structures for 
strain accommodation or prepatterned templates (compliant substrates) [135] 
by rotational alignment. 

5 Conclusions and Outlook 

Molecular dynamics simulations based on empirical potentials provide a suit- 
able tool to study atomic processes that influence macroscopic materials 
properties. The applicability of the technique is demonstrated by calculat- 
ing quantum-dot relaxations and interaction processes with defect creation 
at wafer-bonded interfaces. A brief overview describes the method itself and 
its advantages and limitations, i.e., the macroscopic relevance of the simula- 
tions versus the limited transferability of the potentials. The quality of the 
simulations and the reliability of the results depend on the coupling of the 
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atomistic simulations across length and timescales. Whereas the Lagrange 
formalism is well established for the embedding into suitable environments 
at the continuum level, the approximations of first-principles-based potentials 
need continuous future work. 
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